Testing and tuning symplectic integrators for Hybrid Monte Carlo 

algorithm in lattice QCD 

Tetsuya Takaishi a and Philippe de Forcrand 6,c 
"■Hiroshima University of Economics, Hiroshima 731-0124, Japan 
6 CERN, Physics Department, TH Unit, CH-1211 Geneve 23, Switzerland 
c Institute for Theoretical Physics, ETH Zurich, CH-8093 Zurich, Switzerland 

February 7, 2008 
Abstract 

We examine a new 2nd order integrator recently found by Omelyan et al. The integration error 
of the new integrator measured in the root mean square of the energy difference, (AiJ 2 ) 1 / 2 , is 
about 10 times smaller than that of the standard 2nd order leapfrog (2LF) integrator. As a result, 
the step size of the new integrator can be made about three times larger. Taking into account a 
factor 2 increase in cost, the new integrator is about 50% more efficient than the 2LF integrator. 
Integrating over positions first, then momenta, is slightly more advantageous than the reverse. 
Further parameter tuning is possible. We find that the optimal parameter for the new integrator 
is slightly different from the value obtained by Omelyan et al., and depends on the simulation 
parameters. This integrator could also be advantageous for the Trotter-Suzuki decomposition in 
Quantum Monte Carlo. 

1 Introduction 

The Hybrid Monte Carlo (HMC) algorithm^] is now the established standard for the generation 
of dynamical fermion configurations in lattice Quantum Chromo Dynamics (QCD). The HMC 
algorithm consists of molecular dynamics (MD) trajectories, each followed by a Metropolis test. 
During the MD trajectory, one integrates Hamilton's equations of motion, using an integrator with 
a discrete stepsize At which must satisfy two conditions in order to maintain detailed balance: 
(i) simplecticity (the phase space volume dpdq must be conserved) and (ii) time reversibility. 
The simplest and most widely used integrator wih these properties is the 2nd order leap frog 
(2LF) integrator, which causes 0(Ai 2 ) errors in the total energy or Hamiltonian. These errors are 
eliminated at the Metropolis accept/reject step, which makes the algorithm exact. 

The acceptance at the Metropolis step depends on the magnitude of the error in the total 
energy. In order to reduce the error and thus increase the acceptance one could use higher order 
integrators. Early attempts, however, did not appear to be practical [21 E]. This is because the 
efficiency of higher order integrators depends largely on the system size and these early attempts 
were made on rather small lattices. As the lattice size increases above a certain value V c , the higher 
order integrators should perform better than the low order integrator. This minimum lattice size 
V c depends on the Hamiltonian which we consider and on the choice of integrator. For lattice QCD, 
it turned out that V c becomes very large at small quark masses, so that on currently accessible 
computers the 2LF integrator is the best choice^ for simulations at zero temperature. At finite 
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temperature, higher order integrators could perform better on moderate-size lattices[5]: this is 
because chiral symmetry gets restored, so that small Dirac eigenvalues disappear, which allows for 
stable MD integration using larger stepsizes. 

So far, only the 2LF integrator has been considered in the HMC algorithm of lattice QCD as 
a second order integrator, because of its simplicity and effectiveness. Recently however, Omelyan 
et al. [d] found a new 2nd order integrator which is expected to be better than the 2LF integrator 
although it has twice the computational cost. Here we examine this new 2nd order integrator for 
the HMC algorithm in lattice QCD and measure its efficiency. We also examine the new 4th order 
integrators recommended in [6 . Finally, we try to further tune the new 2nd order integrator. 



2 Symplectic Integrator 

2.1 Recursive construction scheme 

Symplectic integrators are most conveniently described by the Lie algebra formalism [3] E]. Let 
H be a classical Hamiltonian, 

H= l -p 2 + S{q), (f) 

where q = ((ft, (fej • • •) and p = (pi,f>2, • • •) are coordinate variables and conjugate momenta respec- 
tively. Hamilton's equations are expressed as 

/ = {/,#}, (2) 
where / = q or p and {, } stands for the Poisson bracket, 

If we define the linear operator L{H) as 

L(H)f = {f,H}, (4) 

then we can write the formal solution of Hamilton's equations, 

f(t+At) = exp(AtL(H))f(t). (5) 

In general the operator exp(AtL(H)) cannot be expressed exactly in a simple form. Therefore we 
approximate exp(AtL(H)) with an operator correct up to a certain order in At. Let us write L(H) 
as 

L(H) = L(±p 2 )+L(S(q)), (6) 
= T + V, (7) 
where T = L(^p 2 ) and V = L(S(q)). The 2LF integrator is given by decomposing e At ( T + v ) as 

exp(At (T + V)) = exp{^At T) exp(At V") exp^Ai T) + 0(At 3 ). (8) 
We call G 2 (At) the 2LF integrator: 



G 2 {At) = exp(iAt T) exp(Ai V) expf^At T). (9) 
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The integrator G 2 {At) amounts to mapping q and p to new variables as 



q(t + At) \ = f 1 

p(t + At) J I 




) (10) 

= (v 2 (A/).i j. (ii) 



This map is symplectic. This is easy to see, since the three matrices representing the elementary 
substeps are triangular with determinant 1. It is also exactly time-reversible: (^(A^G^— At) = 1. 
An equivalent algorithm is obtained by interchanging T and V in eq. (JSJ . 

Higher order integrators can also be found by decomposing e At ( T + v ) to the desired order. 
Although the decomposition to a higher order is a non-trivial problem with no unique solution, 
there is a simple recursive construction scheme which generates higher order integrators from lower 
order onespfl ITHl In this scheme, the (2k + 2)-th order integrator is given by 

G 2fe+2 (Ai) = G 2k (b 1 At)G2 k {b2At)G 2k {b 1 At), (12) 

where 

6 i = 2 - ' (13) 

2 l/(2fe+l) 

62 = 1 ~ 2 ^ = - 2 -2 WD ' ^ 

Let us call the integrators of ea. l|12|l recursive construction (RC) integrators. These integrators are 
symplectic and constructed in a symmetric way, thus time reversible, i.e. G2fc+2(Ai)G2fc+2(~Ai) = 
1. Note that 62 is negative. The appearance of negative coefficients in higher order integrators is 
inevitable: beyond the 2nd order decomposition there is no decomposition scheme having positive 
coefficients onlv|ll|. However if we include the commutator [V, [T, V]] in the decomposition we may 
circumvent this situation |12l I13| . and integrators with all positive coefficients can be constructed. 
Even so, the inclusion of this commutator requires the calculation of the gradient of the force, which 
increases the computational cost. If the force-gradient calculations are computationally simple for 
the system considered, then it would be worth considering such integrators. For lattice QCD, it is 
unclear that such force-gradient integrators have advantages over the non-force-gradient ones. We 
do not consider such force-gradient integrators here. 



2.2 Minimum norm construction scheme 

Although the recursive construction scheme makes it easy to construct higher order integrators to 
any order, their performance may not be optimal, since the number of force calculations grows 
rapidly with the order of the integrator. More generally, one can decompose e At ( T+v ' as 

exp(Ai [T + V)) = nf exp( Cl At T) exp{d t At V) + 0{At n+1 ), (15) 

where J^- Cj = dj = 1. Moreover, in order to form a time-reversible integrator certain relations 
must be hold. For instance if we take k = 3, the following equations must be satisfied: (i)cj = 
C4, C2 = C3, di = g?3, c?4 = or (m)ci = 0, c 2 = C4, d\ — di, d 2 — d$. For time-reversible integrators, 
the error terms with odd n always vanish|5ll5llllj. Thus time-reversible integrators have a leading 
error term 0(At n+1 ) with n even. 
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The error term 0(At n+1 ) consists of commutators of T and V. For instance, the leading error 
terms of the 2nd and 4th order integrators are respectively 

0(At 3 ) = a[T, [V,T]}+P[V, [VjT]], (16) 

and 

0(At 5 ) = 71 [T, [T, [T, [T, V}}}} + 72 [T, [T, [1/, [T, V}}}} + 73 [V, [T, [T, [T, V}}}} + 74 [V, [V, [T, [T, V}}}}, 

(17) 

where a, (3 and 7$ depend on Ci and d;. 

One strategy to find optimal integrators in the absence of further information about the oper- 
ators T and V is to minimize the norm of the error coefficients. For the case of ea. <|16|l and (|17|l . 
this strategy implies minimizing the following error functions: 

Err 3 = Va 2 +P 2 , (18) 

and 

Err 5 = ^7i+72 +73 +74- (19) 

Omelyan et al. |S] found a class of integrators by following this strategy 1 . Among the new inte- 
grators which they identified, they found several "outstanding" integrators having especially small 
norms of the error coefficients. In this analysis, we consider the new 2nd and 4th order integrators 
which they recommend as outstanding integrators, and which are described as follows. 

2.2.1 2nd order minimum norm (2MN) integrator 

Omelyan et al. [HIE] obtained the following new 2nd order integrator. 

I 2MN {At) = e ^T e ^V e (l-2X)AtT e ^V e XAtT^ (2Q) 

where A takes value A c : 

. 1 (2\/326 + 36) 1 / 3 1 



2 12 {6V326 + 36)V3 

This value of A minimizes \J a{\) 2 + /3(A) 2 , where 



0.1931833275037836. (21) 



„ s 1-6A + 6A 2 
«(A) = T- n , (22) 



12 

- ( 
~24 



/3(A) = ^ (23) 



as can be derived from the expansion of (|20() . 

This integrator requires two force calculations per step. Thus, compared to the 2LF integrator, 
it has twice the computational cost. The norm of the error coefficients Err 3, however, is a factor of 
10 smaller ( Err 2LF / Err 2MN w 10.9 15 ). As we will see later, the error of a 2nd order integrator 
at the end of a Hybrid Monte Carlo trajectory is expected to be proportional to At 2 . Therefore, 
even after taking into account the increased computational cost, we expect that the 2MN integrator 
will perform better than the 2LF integrator, by a factor w \/10.9/2. We will numerically confirm 
this in the next section, and later we will further try to tune the integrator by modifying the error 
function. 



Note that in some cases, the set of polynomial equations defining the optimal decomposition can be solved analytically, 
even beyond the second-order case |14|. 
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2.2.2 4th order minimum norm (4MN) integrator 

At the beginning of the MD integration one can start the integration with either q or p. Usually 
we do not consider this freedom seriously since for the 2nd order integrator the choice of the 
starting variable does not make a significant difference in performance 2 . In general, however, the 
performance could be different depending on the choice of the starting variable. In fact, the optimal 
integrator itself could also be different depending on the starting variable. 

This is precisely what Omelyan et al. found for higher order MN integrators. Let us call velocity 
version 3 the integrator starting by integrating p and position version the integrator starting by 
integrating q. For the optimal 4th order MN integrators with 5 force calculations they found that 
the velocity version has smaller errors than the position version. Actually we have tested both 
integrators and also found that typically the error of the velocity version is a few times smaller 
than that of the position version. In the following numerical tests we use only the velocity version 
of the 4th order MN integrators with 5 force calculations (4MN5FV) which is written as|Hj 

h M N5Fv{At) = e 9Atv e pAtT e XAtv e pAtT e^-^ x+e ^^ v e^- 2 ^ +p ^ AtT 

(24) 

where 

9 = 0.08398315262876693 (25) 
p = 0.2539785108410595 
A = 0.6822365335719091 
fx = -0.03230286765269967. 

Furthermore we have also tested the position version of the 4th order MN integrators with 4 
force calculations (4MN4FP) given by[S| 

h M N4Fp(At) = e P^T e XAtV e 0MT e {l^X)^V e (l-2(9 +P )AtT e (l-2X)^V e 9AtT e XAtV eP AtT^ (2g) 

where 

p = 0.1786178958448091 (27) 

6 = -0.06626458266981843 

A = 0.7123418310626056. 

The velocity version (4MN4FV) is expected to have a similar error to the position version 4MN4FP 
above 6 . Thus we used the 4MN4FP integrator which has one less force evaluation per MD 
trajectory. 



3 Numerical tests of the new integrators 
3.1 Lattice QCD action 

We use the plaquette Wilson gauge and standard Wilson fermion actions with two degenerate 
fermion flavors 16 . The partition function is given by 

Z = J VUdet[M{U)M{U)^]ey^{-S g {U)), (28) 
2 Actually there is a small difference, which we identify later. 

3 One could say momentum version. Here we follow the convention in the literature. 
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with S g (U) the gauge action given by 



s g (u) = ^J2 Tr l 1 - u p( u ^ ( 29 ) 

where U P {U) stands for the plaquette, (3 is the gauge coupling and U is an SU(3) link variable. 
M(U) is the Wilson Dirac operator defined by 

M^U) = S itj +kJ2Ki» - 1)^,m<W - (7p + mUji,j+^ (30) 

where k is the hopping parameter and 7^ are the 7 matrices. The inverse of the hopping parameter 
k is related to the quark mass m and the value k where m is zero is denoted by k c . Asm decreases 
the eigenvalues of the matrix M become small and at the zero quark mass limit the matrix M 
becomes singular. 

Using pseudofermion fields <f> the partition function is re-expressed as 

Z= J VUV^Vcjjexpi-^iM^MiU^y 1 ^- S g (U)}. (31) 

Furthermore introducing momenta p, we obtain 

Z = J VUVtfVcpVp exp[-H(U,(p,^,p)], (32) 

where 

H(C/,^^, p ) = I^;/ + ^t( M ( C/ ) M ( C/ )t)-^ + 5 s ( C ;) (33) 
is the Hamiltonian we consider in our numerical tests. 



3.2 Hybrid Monte Carlo algorithm 

The HMC algorithm combines MD and Metropolis accept/reject stepspQ to form a Markov chain. 
Starting from an "old" configuration {[/}, a "candidate" configuration {[/'} is obtained by (i) 
drawing momenta {p} and pseudofermion fields {0} from Gaussian distributions exp(— \p 2 ) and 
exp(—<p^ (M(U)M(U) ji )~ 1 (j)) respectively; (ii) integrating Hamilton's equations of motion cq.© 
with a discrete stepsize integrator. In order to maintain detailed balance the integrator in the 
MD step must satisfy two conditions: simplecticity and time reversibility. Let TM_o(At) be an 
elementary MD step with a discrete step size At. Tj^n^At) evolves (p, (?) to (p',q'): 

T MD (At):(p,q)^(p',q>). (34) 

The time reversibility condition means the following is satisfied. 

T M D(-At):(p',q')^(p,q). (35) 

The simplecticity is the condition that the phase space volume is conserved, dpdq = dp' dq' . The 
symplectic integrators in Sec. 2 satisfy the above two conditions. 

This elementary MD step is performed repeatedly to integrate the equations up to a certain 
"length" of the MD trajectory. Here the trajectory length is set to unity. In general the integrator 
can not solve the Hamilton's equations of motion exactly and thus the energy is not conserved. Let 
AH be the change in energy caused by the integrator at the end of the trajectory. This error is 
corrected by the Metropolis test, which accepts the candidate configuration {U'} at the end of the 
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trajectory with probability p = min(exp(— AH), 1). If the candidate configuration is rejected, the 
old configuration {U} is included again in the Markov chain. 

Thus, a trade-off must be achieved between two conflicting goals: good energy conservation, 
because it ensures high acceptance probability in the Metropolis step, and low computer cost. This 
is why the choice of integrator plays a crucial role. 

The fermionic part of the force dH/dU is given by 

dH f^~ = -^(MJtft)- 1 9 ^^ (MMt)-V- (36) 

Since the matrix M is singular at the quark mass m — 0, as m decreases the fermionic force 
diverges. Thus HMC simulations in lattice QCD become unstable at small quark masses, unless 
the stepsize is reduced in proportion to m. 



3.3 Error of Hamiltonian and Acceptance 

Here we summarize the expected behavior for the error of the Hamiltonian and the acceptance 
of the HMC algorithm. The n-th order integrator causes 0(At n+1 ) integration errors for q and 
p. However the error in the Hamiltonian at the end of a unit-time trajectory 4 is 0(At n ). Thus 
AH ~ At n . Furthermore from Creutz's equality (exp(Aff)) = 1 JSj we expect 

(AH) ~ (AH 2 ) cx VAt 2n , (37) 

where V is the volume of the system. Thus the root mean square of the error of the Hamiltonian 
at small At is expected to be 

(AH 2 ) 1 ' 2 « C n V 1/2 At n , (38) 

where C n is a Hamiltonian-dependent coefficient. Using (AH 2 ) 1 ' 2 , the acceptance of the HMC 
algorithm for large volumes is given by |17| 



(P 



<!<<■/ 



'8 

For small (AH 2 ) 1 ' 2 , one may use the approximate formula 

1 



erfc^Atf 2 ) 1 / 2 ). (39) 



exp(- — (AI/ 2 ) 1 ^), (40) 



which is applicable for P acc > 20% /JJ. The performance of integrators can be measured by the 
[inverse of the] work per accepted trajectory, i.e. by the product of the acceptance and step size: 
P acc x At. The best performance of integrators is obtained at the step size which maximize P acc X At. 
Using eqs.lJSSJ and we obtain the optimal acceptance which maximizes P acc x At as 0] 

P opt = exp f-ij , (41) 

which does not depend on the details of the Hamiltonian but only on the order of the integrator. 
This result indicates that the optimal acceptance for any 2nd order integrator is about 61% which 
is consistent with the numerical results of 60 ~ 70% 4 . Ea. (|41() also indicates that the optimal 
acceptance increases with the order of the integrator: 78% for 4th order and 85% for 6th order. 



4 AH does not increase linearly with trajectory length. It increases linearly up to a certain characteristic length l c 
provided that At is not too large, then saturates. Thus the accumulated error in the Hamiltonian is expected to be 

la 
At 



AH ~ 0(At n+1 ) x i- ~ Q(At n )\n\. 
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3.4 Performance of 2nd order MN (2MN) integrator 

Here we compare the efficiency of the 2MN integrator with that of the 2LF integrator. For 2nd 
order integrators, from eq. P%| (AH 2 ) 1/2 at small At is expected to be C 2 V 1 / 2 At 2 . We measure 
the coefficient C 2 for both integrators at small enough At, and by comparing the coefficients we 
obtain the performance of the 2MN integrator relative to the 2LF integrator. 

Figure^shows (Ai? 2 ) 1 / 2 as a function of step size At at P = 5.00 and k = 0.160 on 4 4 lattices. 
We see that (AT? 2 ) 1 / 2 is proportional to At 2 as expected and the error of the 2MN integrator is 
about 10 times smaller than that of the 2LF integrator at any At until instabilities show up. 

Figure |21 shows the ratio C 2 lf/C 2 mn as a function of k. The coefficients C 2 lf and C 2 mn are 
extracted by using eci. (|38H for the 2nd order with simulations at a small value of At. As seen in 
the figure, C 2 lf /C 2 mn is about 10, which means that the error of the 2MN integrator is about 10 
times smaller than that of the 2LF integrator. This is consistent with the theoretical expectation of 
Omelyan et at. If we take C 2 lf /C 2 mn ~ 10, this means that the step size of the 2MN integrator 
can be increased by a factor 3 (~ VlO) over that of the 2LF integrator, as long as the error still 
behaves as At 2 . Since the 2MN i ntegrator has two force calculations per elementary step, the 
efficiency should be measured by y C 2 lf/C 2 mn /2, which is about 1.5. Thus it is concluded that 
the 2MN integrator is about 50% faster than the 2LF integrator. 

Of course, this assessment rests on the assumption that the step size can indeed be increased 
without running into instabilities, so that the limiting factor in the step size comes from the error 
accumulation. Note that the 2MN integrator appears no worse, or perhaps slightly better, than 
the 2LF in terms of instabilities: departure from the quadratic behaviour (AiJ 2 ) 1 / 2 cx At 2 starts 
at similar values of (Aff 2 ) 1 / 2 in Fig.^ and appears more gradual. 

3.5 Comparison of 2nd and 4th order MN integrators 

The efficiency of higher order integrators should be measured against lower order ones. From the 
above analysis we know that the 2MN integrator is more efficient than the 2LF integrator. Therefore 
we compare the 4MN integrator with the 2MN integrator. Assuming eq. I|38l) the comparison could 
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Figure 2: C2LF/C2MN as a function of k. Simulations at /3 = 5.00(5.60) are performed on 4 4 (8 4 ) 
lattices. 



be done by following the analysis of @|. However we found a problem with the 4MN integrator. 
Namely the error of the Hamiltonian (AH 2 ) 1 / 2 is not simply described by ea. H38l) but is dominated 
by higher order terms in At already at small (AH 2 ) 1 / 2 . Figure [3] shows (AH 2 ) 1 / 2 on 8 4 lattices as a 
function of At. As seen in the figure (left), at a fixed step size the error of the 4MN5FV integrator is 
about 1000 times smaller than the previously known 4th order integrator (4RC), which is consistent 
with the theoretical expectation^. The expected behavior of (AH 2 ) 1 / 2 w C±V 1 / 2 At 4: , however, is 
seen only at small (AH 2 ) 1 / 2 . We are only interested in the region of 0.1 < (AH 2 ) 1 / 2 < 1 which 
corresponds to an acceptance of 60% ~ 95% 5 . In this region, (AH 2 ) 1 / 2 of the 4MN5FV integrator 
is dominated by higher order terms in At, which results in that (AH 2 ) 1 / 2 grows rapidly with At. 
This observation makes the 4MN integrator unattractive on practical lattice sizes. 

Although the 4MN4FP integrator seems to be more stable than the 4MN5FV integrator, it also 
shows the deviation from the At 4 line at small (AH 2 ) 1 / 2 (See Fig. 3). Thus compared to the 2MN 
integrator, the 4MN integrators tested here seem unattractive. As the quark mass m decreases, 
it is often seen that the integrator becomes unstable |2()|. because the force increases as 1/m. In 
lattice QCD calculations the parameter region of small quark masses is the physically interesting 
one. In this region the 4MN integrators may easily show instability, which limits their applicability. 

At finite temperature, however, the coefficients C n behave differently from those at zero tem- 
perature. Typically we expect < Cj =0 . Therefore at finite temperature we may be lead to 
a different conclusion and this must be studied numerically. A numerical test showed that at finite 
temperature the 4RC integrator performs better than the 2LF integrator on lattices larger than a 
minimum sizej^]. We have made the same test for the 2MN and 4MN5FV integrators on an 18 3 x 4 
lattice at j3 = 5.75 and n — 0.1525. For the 2MN integrator the acceptance is measured to be about 
0.6 at At = 0.1 and for the 4MN5FV integrator the acceptance is about 0.8 at At = 0.37. These 
values of the acceptance are close enough to the optimal acceptance given by ea. lj4*T)l . The gain of 
the 4MN5FV integrator over the 2MN one is calculated by 

17= , (42) 

5 Figure 1 of g] 
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At At 



Figure 3: (AH 2 ) 1 / 2 as a function of At. Simulations are performed on 8 4 lattices. 4RC stands for the 
4th order integrator obtained by ea. Q12|l . The lines proportional to At n are drawn to guide the eye. 



where K42 is the relative cost factor and K42 = 2.5 for the 2MN and 4MN5FV integrators. Sub- 
stituting the measured values into G, G is calculated to be ss 2, which shows that the 4MN5FV 
integrator is more effective. Thus at finite temperature there is room to use a 4MN integrator 
depending on the simulation parameters. 

4 Tuning the 2MN integrator 

The strategy to minimize eas. <|18[) and (|19fl is based on the assumption that the errors coming from 
the two commutators [T, [V, T]] and [V, [V,T]] are equally dominant. In general this simplifying 
assumption does not hold. Here we try to minimize a more general form of the error function. Let 
us assume the following form of (AH 2 ) 1 / 2 , 

(AH 2 ) 1 / 2 « ^a(\) 2 f{At) 2 +(3(\) 2 g(At) 2 , (43) 
= V«(A) 2 / 2 + /?(A)VAt 2 , (44) 

where a(A) and /3(A) are given by eas. (|22l) and and f 2 and g 2 are unknown parameters, to be 
determined from numerical simulations. In general, by performing simulations at two values of A 
one can determine f 2 and g 2 numerically. The determination can be made easier by noticing that 
a(Ai) and (3(\2) are zero at Ai = (1 — l/\/3)/2 and A2 = 1/6 respectively. By simulating at Ai and 
A2 we immediately obtain f 2 and g 2 . 

Figure^lshows (AH 2 ) 1 / 2 at (3 — 5.00 and k = 0.160 as a function of A. We see that the optimal 
A at the minimum of (AH 2 ) 1 / 2 is slightly different from A c , the value of eci. (|21J) . Moreover the 
optimal A is different between the velocity and position version integrators. Here ea. (|20J) is the 
position version integrator. The velocity version is obtained by interchanging T and V in ea. (|20|l . 
The lines in the figure are given by ea . Q44[l . with f 2 and g 2 determined by simulations of the 
position version of 2MN integrator at Ai and A2. To draw the dashed line for the velocity version, 
we simply interchange the values f 2 and g 2 . Both lines describe the numerical results very well, 
down to A = 0. Note that the velocity version of 2MN integrator becomes the position version of 
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Figure 4: (AH 2 ) 1 / 2 as a function of A. The right figure is a zoom of the left. Simulations are 
performed at /3 = 5.00 and n = 0.160 on 4 4 lattices with At = 0.05. The lines are determined from 
simulations of the position version integrator at Ai and A2. The position version has a small advantage 
over the velocity version, since it gives a slightly reduced minimum RMS error (right). 



2LF integrator at A = 0, and vice versa. The position version of the 2MN integrator gives a slightly 
smaller minimum. At A = the velocity version of 2MN integrator has a smaller error than the 
position version, which means that the position version of 2LF integrator has a smaller error than 
the velocity version. This was already observed in |21j . Since the position version also leads to one 
less force evaluation by the end of a trajectory we definitely recommend using the position version 
(for the 2MN and the 2LF integrators both): it requires less work and gives a higher acceptance. 

The quality of our fit justifies a posteriori the ansatz made for the magnitude of the error ea . (|44|1 . 
Indeed, to leading order At 4 , the error (AH 2 ) should be of the form 

(AH 2 ) = (a{X) 2 (F 2 ) + a{X)(3(X)(FG + GF) + /3(A) 2 (G 2 ))At 4 , (45) 

where F = [T, [V, T]] and G = [V,[V,T]] from eq.©. We find that the crossterm (FG) is in 
our case "one order of magnitude smaller" than (F 2 ) and (G 2 ), indicating that the two operators 
are almost uncorrelated in our system. While this finding may not be true in general, it provides 
support for the minimum norm strategy of Omelyan et al. at least in the context of lattice QCD. 

Figure shows \f\At 2 , \g\At 2 and |/|/|g| as a function of 1/k at — 5.00 on 4 4 lattices 
(At = 0.05). As one approaches k c , both |/| and \g\ increase. On the other hand the ratio 
decreases. This is as expected since g comes from the error term [V, [V,T]] involving two factors of 
the fermion potential, versus one for / and [T, [V, T]]. Although |/|/|<?| seems to approach one at n c , 
there is a possibility that it further goes down to zero and its limit must be carefully investigated. 
Note that when |/|/|<7| = 1, the optimal A becomes A c . 

Figure El shows the optimal A as a function of We see that the optimal A is different from 
A c and slightly larger. 

5 Conclusions 

We have tested the new 2nd and 4th order integrators obtained by minimizing the norm of the error 
coefficients. We find that the 2MN integrator performs better than the conventional 2LF integrator, 
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Figure 5: |/ |At 2 , |g|At 2 and |/|/|g| as a function of 1/k. Simulations are performed at (3 
4 4 lattices with At = 0.05. The dashed line indicates k c = 0.187["H?]. 
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Figure 6: Optimal A as a function of 1/k. Simulations are performed at (3 = 5.00 on 4 4 lattices with 
At = 0.05. The dashed line indicates k c = 0.187. 
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by about 50%. Therefore we recommend to use the 2MN integrator in HMC simulations. Although 
in our tests we used the standard Wilson fermion action, the 2MN integrator can be used for any 
actions, e.g, KS fermions, improved actions, polynomial actions for odd flavors |2*2l 123) . Moreover 
we may combine the 2MN integrator with other acceleration techniques such as multiple time step 
integration [5], multiple pseudo- fermions 24: and preconditioned actions [25]. 

The same 2MN integrator can also be used in the Trotter-Suzuki decomposition of the partition 
function: exp(— f3H) = (exp(— At H)) , where N — f3/At, in Quantum Monte Carlo simulations, 
when a formulation in continuous imaginary time |25| is not practical. 

Although at first sight, one can equivalently start by integrating over positions or velocities, we 
observe that integrating over positions first gives a slightly higher acceptance |2*T] . with one less 
force evaluation at the end of a trajectory. 
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